#include "Coils.h"
#include "GreenFunction.h"
#include "Log.h"
#include "Parameters.h"
#include "Grid.h"
#include "DynamicArray.h"
#include "Psi.h"

#include <cmath>
#include <math.h>
#include <iostream>
#include <stdio.h>

extern Log* plog;
extern Parameters* ppmts;

Psi::Psi(Grid* grid, Coils* coils) {
     num_r = grid->getNumberR();
     num_z = grid->getNumberZ();
     right_r = grid->getRightR();
     left_r = grid->getLeftR();
     up_z = grid->getUpZ();
     down_z = grid->getDownZ();
     h_r = grid->getHR();
     h_z = grid->getHZ();
     psi_g = ppmts->getPsiG();
   
     num_coil = coils->getNumber();
     its_coils = coils;
     delta_s = grid->getDelta();

     //initalizing plasma area indexes
    left_r_g = (ppmts->getPlasmaLeftR() - left_r) / h_r;
    right_r_g = (ppmts->getPlasmaRightR() - left_r) / h_r;
    up_z_g = (ppmts->getPlasmaUpZ() - down_z) / h_z;
    down_z_g = (ppmts->getPlasmaDownZ() - down_z) / h_z;

    m = up_z_g - down_z_g;
    n = right_r_g - left_r_g;
     psist = AllocateDynamicArray<float>(num_r, num_z);
     psipl_1 = AllocateDynamicArray<float>(num_r, num_z);
     psipl_m = AllocateDynamicArray<float>(num_r, num_z);
     psipl_tmp = AllocateDynamicArray<float>(num_r, num_z);

     psipl_s = AllocateDynamicArray<float>(num_r, num_z);
     psipl_rs = AllocateDynamicArray<float>(num_r, num_z);
     psipl_on_g_s = AllocateDynamicArray<float>(num_r, num_z);
     psitotal_s = AllocateDynamicArray<float>(num_r, num_z);
     psitotal_1 = AllocateDynamicArray<float>(num_r, num_z);
     psitotal_m = AllocateDynamicArray<float>(num_r, num_z);
     psitotal_tmp = AllocateDynamicArray<float>(num_r, num_z);

     psipl_on_g = AllocateDynamicArray<float>(num_r, num_z);

};

Psi::~Psi(){
    FreeDynamicArray<float>(psist);
    FreeDynamicArray<float>(psipl_1);
    FreeDynamicArray<float>(psipl_m);
    FreeDynamicArray<float>(psipl_tmp);


    FreeDynamicArray<float>(psitotal_1);
    FreeDynamicArray<float>(psitotal_m);
    FreeDynamicArray<float>(psitotal_tmp);
    FreeDynamicArray<float>(psipl_on_g);

    FreeDynamicArray<float>(psipl_s);
    FreeDynamicArray<float>(psipl_rs);
    FreeDynamicArray<float>(psitotal_s);
    FreeDynamicArray<float>(psipl_on_g_s);


};

float Psi::st_point(float r, float z) {
    float res = 0;

    for (int i = 0; i < num_coil; i++) {
        res = res + G(r, z, its_coils->getCoilR(i), its_coils->getCoilZ(i)) * its_coils->getCoilI(i);

    }
    if (res != res) {
        plog->writeErr("psi_st");
    }
    return res;
}

float** Psi::st_grid() {

    float r, z = 0;
 
    for (int j = 0; j < num_z; j++) {
        z = down_z + j*h_z;
        for (int i = 0; i < num_r; i++) {
            psist[i][j] = 0;
            r = left_r + i*h_r;
           
            psist[i][j] = st_point(r, z);

        }
    }
    return psist;

};

float Psi::pl_point(float r, float z, float** j_phi) {
    float res = 0;
    float ri, zi = 0;
    for (int j = 0; j < num_z; j++) {
        zi = down_z + j*h_z;
        for (int i = 0; i < num_r; i++) {
           
            ri = left_r + i*h_r;
          
            res = res + G(r,z, ri, zi)*j_phi[i][j];

        }
    }
    res = res*delta_s;
    return res;


};

float** Psi::pl_grid(float** j_phi, float** psipl) {

    float r, z = 0;

    for (int j = 0; j < num_z; j++) {
        z = down_z + j*h_z;
        for (int i = 0; i < num_r; i++) {
          
            r = left_r + i*h_r;

            psipl[i][j] = pl_point(r, z, j_phi);

        }
    }
    return psipl;
};


float Psi::total_point(float psi_pl, float psi_st) {
    return psi_pl+psi_st;

};

float** Psi::total_grid(float** psi_pl, float** psi_st, float** psisum) {

    for (int j = 0; j < num_z; j++) {

        for (int i = 0; i < num_r; i++) {


            psisum[i][j] = psi_pl[i][j] + psi_st[i][j];

        }
    }
    return psisum;

}

void Psi::bound(float** psi) {
    FILE* out_file;
    out_file = fopen("bound.txt", "wb+");
    float r, z = 0;
    for (int j = 0; j < num_z; j++) {
        z = down_z + j*h_z;

        for (int i = 0; i < num_r; i++) {
           r = left_r + i*h_r;
           if ( abs(psi[i][j] - psi_g) < ppmts->getPrecision() ) {
               fprintf(out_file, "%f %f %f\n", r, z, psi[i][j]);
           }


        }
    }
    fclose(out_file);

};


float** Psi::get_psipl_on_g(float** jphi) {
    float res = 0;

    float r, z = 0;
    float ri, zi = 0;
    for (int j = down_z_g; j <= up_z_g; j++) {
       

       for (int i = left_r_g; i <= right_r_g; i++) {
          if ( i ==  left_r_g || i == right_r_g || j == down_z_g || j == up_z_g ) {
                r = left_r + i * h_r;
                z = down_z + j * h_z;
                /////////////////// integral calculation
                res = 0;
                ri, zi = 0;
                for (int m = down_z_g; m <= up_z_g; m++) {
                    zi = down_z + m * h_z;

                    for(int n = left_r_g; n <= right_r_g; n++){
                        ri = left_r + n * h_r;
                        res = res + G(r, z, ri, zi)*jphi[n][m];

                    }

                }
                res = res * delta_s;
                psipl_on_g[i][j] = res;                       
          }
       }
        
    }
    return psipl_on_g;

};

float Psi::H(float psi_) {
    if ( psi_ < 0) {return 0;}
    else {return 1;}

};

float Psi::F(int i, int j) {

};